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A New Multidimensional Relativistic Hydrodynamics code based 
on Semidiscrete Central and WENO schemes. 

Tanvir Rahman Robert. B. Moore^ 

ABSTRACT 

We have proposed a new High Resolution Shock Capturing (HRSC) scheme 
for Special Relativistic Hydrodynamics (SRHD) based on the semidiscrete cen¬ 
tral Godunov-type schemes and a modified Weighted Essentially Non-oscillatory 
(WENO) data reconstruction algorithm. This is the first application of the 
semidiscrete central schemes with high order WENO data reconstruction to the 
SRHD equations. This method does not use a Riemann solver for flux computa¬ 
tions and a number of one and two dimensional benchmark tests show that the 
algorithm is robust and comparable in accuracy to other SRHD codes. 


1. Introduction 

Gas flows at (ultra) relativistic speeds are an integral component of many astrophysi- 
cal phenomena. Some significant examples of these include accretion disks around compact 
objects, Gamma ray bursts, collapse of (super)massive stars, pulsar wind nebula/supernova 
remnant interactions. Active Galactic Nuclei (AGN), X-ray binaries, superluminal jets and 
many others. As it stands, there is a massive amount of astronomical data from many of 
these sources that to date, have only been partially understood and which hold clues to 
many problems related to the phenomena mentioned above. Gorrelating these data with 
the theory in many of these phenomena would require doing simulations that would include 
modelling (ultra)relativistic gas/fluid flows. Such flows are described by non-linear hyper¬ 
bolic conservation laws that have shock waves as possible solutions. Shock waves are difficult 
to approximate numerically using standard finite difference techniques as they usually lead 
to spurious oscillations. High Resolution Shock Gapturing (HRSG) schemes are a class of 
numerical methods devoted specifically to approximating hyperbolic conservation laws. Over 
the past fifteen years HRSG research have progressed immensely and they have been applied 
to a wide variety of problems involving classical and relativistic hydrodynamics. For a ped¬ 
agogical introduction to some modern HRSG schemes and their applications, see Levequque 
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R. (2002); Toro (1999); Pen U. L. (2002). To date, the most commonly used HRSC 
schemes in relativistic hydrodynamics have been ones that have not incorporated many of 
the recent advances in HRSC research. In this work, we will apply a new scheme to the 
SRHD equations that is simpler than most of the previous approaches and incorporates a 
number of recent advances from HRSC research. In Rahman & Moore (2005), this scheme 
was applied to the multidimensional classical hydrodynamics code. The work presented here 
should be considered an extension of that work. Before describing the motivations behind 
the algorithm used in this work, we provide a brief review of numerical methods used for the 
SRHD equations and that of HRSC schemes. 

Wilson J. R. (1972), Wilson et al. (1979), Hawley et al. (1984), and Centrella & Wilson 
(1984) were the hrst to apply numerical methods to approximate the SRHD equations. 
They discretized the SRHD equations using an explicit hnite difference scheme with artihcial 
viscosity and monotonic transport. Since their work, this technique has been used to study 
a number of astrophysical phenomena including stellar collapses, accretion disks, cosmology 
etc. Most SRHD schemes up to the late eighties used the artihcial viscosity technique of 
Wilson et al. (1979) to handle shock waves. However, a major breakthrough occurred when 
HRSC schemes were applied to the SRHD equations (Marti et al. (1991); Marquina et. al 
(1992); Eulderink, F (1993, 1995)). One of the most signihcant of these new schemes was 
the relativistic Piecewise Parabolic Method (PPM) by Marti & Muller (1996) (MM) which 
was based on the PPM reconstruction scheme of Colella & Woodward (1985). For a recent 
review of modern SRHD schemes, we refer the reader to Marti & Muller (2003). 

The use of HRSC schemes for SRHD equations essentially shifted the paradigm for 
solving these equations from artihcial viscosity based methods to this approach. Modern 
HRSC schemes are based on two conventional approaches for solving hyperbolic conservation 
laws. These are the so-called upwind (Harten A. (1983); Van Leer B. (1979); Roe P. L., 
(1981); Toro (1999)) and central methods (Lax & Friedrichs (1971); Lax P. (1954); Nessyahu 
& Tadmore (1990)). A common aspect of most modern HRSC-SRHD schemes is that they 
are solved using upwind schemes that use so-called Riemann solvers to approximate the 
solution. A drawback of using Riemann solvers is that they are computationally expensive 
and difficult to implement because they need computations of eigenvalues and eigenvectors 
of the hux matrix and uses hux splittings etc., to compute the numerical huxes. On the other 
hand, central schemes do not use Riemann solvers and are simpler and easier to implement. 
Until recently, central schemes were not considered for wide ranging applications because 
a trade-off for their simplicity was loss of accuracy. This is because central schemes were 
generally more dissipative than their upwind counterparts. However, recent progress in 
HRSC research has addressed this issue and some newer formulations of central schemes 
have been shown to be comparable performance-wise to the upwind approach. The work 
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presented here is based on some of these newer central schemes. Before discussing them, we 
provide below a brief overview of central schemes and discuss the most important aspect of 
HRSC schemes known as non-oscillatory data reconstructions methods. 

Modern central schemes are based on the first and second order shock capturing algo¬ 
rithms developed by Lax-Friedrichs (LxF) (Lax & Friedrichs (1971); Lax P. (1954)), and 
Nessyahu & Tadmore (1990) (NT), respectively. Since NT, there have been a number of 
extensions of central HRSC schemes to higher orders and multidimensions. Some of these 
recent extensions are summarized below. The NT scheme was extended to 3’’'^ order by 
Liu & Tadmor (1998), and to multidimensions by Jiang & Tadmore (1998). These were fol¬ 
lowed by a new formulation of central HRSC schemes known as semidiscrete central schemes. 
The semidiscrete schemes were designed to address the dissipation issue mentioned above 
and were proposed in Kurganov & Tadmore (2000) (KT). KT showed that semidiscrete 
schemes, which are formulated on non-staggered grids and use more accurate information 
of local speeds of propagation, are less dissipative than their staggered grid counterparts. 
Hence, these schemes retained the simplicity of the central approach and were comparable 
in accuracy to other Riemann solver based upwind approaches. The inception of semidis¬ 
crete central schemes precipitated a great deal research on extending this formalism which 
included higher order formulations in multidimensions, genuinely multidimensional formu¬ 
lations and unstructured grid formulations among others(Liu & Tadmore (1998); Jiang & 
Tadmore (1998); Levy et al. (2000, 2002); Kurganov &Petrova (2001); Kurganov & Noelle 
(2001); Kurganov & Levy (2000)). 

Besides the method used to advance the solution, another aspect of HRSC schemes that 
is equally important for any algorithm is the so called non-oscillatory data interpolation 
method. Data interpolation is used in all HRSC schemes. Non-oscillatory data reconstruc¬ 
tion is the piecewise continuous polynomial interpolation of the data (which may contain 
discontinuities) over the computational domain. Most of the HRSC schemes mentioned 
above use the so called Essentially Non-Oscillatory (ENO) (Harten et al. (1987)) data 
reconstruction algorithm for data interpolation. This methods works by interpolating the 
data using the smoothest stencil from a number of choices. A modern extension of ENO 
schemes is the Weighted Essentially Non-Oscillatory data reconstruction (WENO) algorithm 
(Liu et al. (1994); Jiang G-S, Shu C-W (1996)) that has a number of advantages over its 
predecessor. In order to take advantage of the developments in central approaches and high 
order data reconstruction, another class of HRSC schemes were developed as an extension of 
the NT scheme by Levy et al. (Levy et al. (2000, 2002)) called the Central Weighted Non- 
Oscillatory (CWENO) schemes. Recently, Kurganov & Levy (2000) (KL) have combined 
the semidiscrete central schemes with the WENO reconstruction method and proposed yet 
another, better HRSC scheme. The work presented here is based on the KL HRSC scheme 
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which we will describe below. 

Despite the many advances mentioned above, there have been relatively few attempts 
at applying central-type methods to computational astrophysics. However this is rapidly 
changing. Among the most significant attempts to apply central schemes for computational 
astrophysics, the following are noteworthy. Del Zanna & Bucciantini (2002), have developed 
an algorithm that is a variation of the conventional finite volume central type approach and 
applied it to multidimensional SRHD and SRMHD equations. They have also used these 
codes to study Pulsar Wind Nebula (PWN) interactions with the interstellar medium Del 
Zanna & Bucciantini (2004) and supernova remnants Del Zanna & Bucciantini (2003). 
Anninos & Fragile (2003,2004) have developed a new central type scheme for relativistic 
hydrodynamics in fixed, curved spacetime and have applied their scheme to study accretion 
disks around kerr black holes. Lucas-Serano et ah (2004) have assessed the applicability of 
central schemes to SRHD equations using the PPM data reconstruction scheme Colella & 
Woodward (1985). Their results were the first to demonstrate that central type schemes of 
KT can be used for SRHD. Recently, Shibata & Font (2005) have successfully tested the suit¬ 
ability of central type schemes for general relativistic simulations. We turn our attention now 
to applications of the WENO data reconstruction methods in computational astrophysics. 
Balsara (Balsara D. (2001)) have used the WENO methods extensively in magnetohydro¬ 
dynamics. Feng et ah (2004), have recently developed a cosmological hydrodynamics code 
using the WENO reconstruction scheme. Zhang & MacFadyen (2005), have recently used 
the WENO scheme in their adaptive SRHD code. The work mentioned above point to a 
promising future for the application of central and WENO schemes in astrophysical hydro¬ 
dynamics. Indeed, the work presented here was done concurrently with many of the recent 
work mentioned above. However, the algorithm presented here is the first to combine the 
WENO and the semidiscrete central approach for applications in SRHD. We describe below 
this algorithm and its novel aspects. 

In a previous paper, Rahman & Moore (2005), we have described in detail the in¬ 
centive behind our algorithm for multidimensional classical hydrodynamics. Many of the 
same reasons apply for the SRHD equations as well. The algorithm is based on combin¬ 
ing the semidiscrete approach of Kurganov & Levy (2000) with high order WENO data 
reconstruction methods. To ensure robustness we have added the flattening, steepening and 
monotonicity preserving techniques of the PPM reconstruction scheme by Colella & Wood¬ 
ward (1985) to the data reconstruction scheme. Essentially, the simplicity of central schemes, 
the accuracy WENO reconstructions and the robustness of the piecewise parabolic method 
have been combined to propose a new robust central scheme. Building on the success of this 
algorithm for classical hydrodynamics, it has been applied to the SRHD equations. This 
work is original in several respects. It is the first attempt at solving the SRHD equations 
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using the semidiscrete central scheme with WENO data reconstructions. The algorithm 
has been tested using both and 4*^' order data reconstruction techniques. It is also a 
dimensionally unsplit algorithm and to our knowledge is the only unsplit multidimensional 
algorithm in SRHD. The work presented here can be considered as laying the foundation for 
the development of a multi-purpose, multi-dimensional SRHD code that could be used to 
study a wide variety of astrophysical phenomena. 

The paper is organized as follows. Sec. 2 presents the SRHD equations. Sec. 3 describes 
the new algorithm for numerically approximating the SRHD equations. In Sec. 4, we present 
the results of a number of one dimensional benchmark tests. In Sec. 5, the two dimensional 
tests are presented. Finally we conclude this paper in Sec. 6 . 


2. Relativistic Hydrodynamic Equations 

As a system of hyperbolic conservation laws, the equations of special relativistic 
dynamics are given by 

(9U(w) (9f*(w) 

--^ = 0 . 

dt dx^ 

where the indices run from 1 to 3. In the above equations U and f* are given by 


U(w) = (D ,S‘ ,t) , 

(2) 

r(w) = {Dv‘ ,SV+p5« ,S‘-Dv‘) . 

(3) 

Here 6^^ is the Kronecker delta, n* is the three velocity and p is the pressure, 
between the conserved variable and the primitive variable is given by 

The relationship 

D = pW , 

(4) 

= phW\^ , 

(5) 

T = phlV"^ — p — D , 

( 6 ) 

where h is the specific enthalpy; h = 1 + e + p/p, and W is the Lorentz factor satisfying 
W = vP = 1/Vl — with = v^Vi. The 3-velocity components are obtained from the 

spatial components of the 4-velocity as n* = /u^. The equations above 

ideal gas Eqnation of State (EOS) given by. 

are closed by the 

P = (7 - l)pe , 

(7) 


hydro- 

( 1 ) 


where 7 is the adiabatic index of the ideal gas and e the energy density. The difference be¬ 
tween the conserved and the primitive variables in these equations require special treatment 
in the algorithm which will be discussed below. 
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3. The Semidiscrete Central WENO Algorithm for SRHD 


This section presents the new algorithm for solving the SRHD equations implemented in 
this work. This algorithm is very similar to that presented in a previous paper (?). Besides 
one major sub-step, that of recovering the primitive variables from the conserved variables, 
it is in fact identical. However, for the sake of completeness, a very brief description of it is 
provided. 


As mentioned before, this new SRHD algorithm is based on the combination of the 
Semidiscrete Central scheme of Kurganov & Levy (2000) (KL), and the steepening, flattening 
and monotonicity preserving algorithm of MM. For a review and details of the derivation 
of KL, see Kurganov & Levy (2000). The reader is also referred to MM for details of the 
steepening, flattening and monotonicity preservation algorithms. In KL, the solution of a 
scalar hyperbolic conservation law in one dimension. 


Ut + 


dfju) 

dx 


= 0 


u{x, t = t^) = u^{x) 


( 8 ) 


is given by 




-ruAt) = -^^ 

dt ^ Ax 


where the flux H-,i is. 


:=-^ 


r 




uT^i, Oj.,.! are given by 


“i+i = q‘(a:j+i/2) , 


(9) 


( 10 ) 


( 11 ) 


“i+J :=m£“{/>(^(“,+i/2)).P(^(“t+i/2))} . (12) 

where {P”(a;)} is a non-oscillatory piecewise polynomial reconstruction of cell averaged data 
at time t = t”, {Uj{x)}. In this paper, both and 4*^ order WENO reconstruction methods 
of Kurganov & Levy (2000); Levy et al. (2000) are used. In Eq. 10, is the speed 
of propagation of u at the interface of a cell that is determined from the spectral radius of 
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the Jacobian of the flux /. For a system of equations in multidimensions, the scheme given 
above can be extended as follows. A dimensionally unsplit approach can be used to extend 
Eq. 9 to multidimensions. In this approach, Eq. 9 is modihed by adding another term 
corresponding to the flux differences in the y-direction. For a system of equations, the scalar 
/ is replaced by a vector F and is replaced by the Jacobian of F. 

For all the computations presented here, the piecewise polynomial interpolations have 
been done on the primitive variables from which the conserved variables are computed using 
Eqns. 4-6. In each step of our integration routine, once the conserved variables are updated, 
the primitive variables are immediately recovered. The method used to do this is by solving 
of a non-linear equation using the Newton-Raphson root Ending method (see Marti & Miiller 
(2003) for details). It is this aspect of numerical SRHD that makes them computationally 
more expensive than the Euler equations of gas dynamics. We can now summarize our 
algorithm in the following few steps: 

step 1; Given initial data u" (primitive variables), use an nth. (for n < 4 ) order 
WENO reeonstruetion algorithm to construet a piecewise polynomial interpolation of each 
variable. 

step 2; Apply the steepening (only to the density p), flattening and monotonicity pre¬ 
serving algorithms to the interface values of the primitive variables. 

step 3: Update Ttf (the conserved variables) to u!)^^ using Egns. 9 and 10 above. 

step 4: Recover the primitive variables from the updated conserved variables. 

step 5: Repeat steps 1 to step 4 

The time integration of Eq. 9 is done using a high order total variation diminishing 
(TVD) Runge-Kutta scheme Shu & Osher (1988), which combines the first order forward 
Euler method with predictor-corrector steps. This is, 

jjii) = un^ ATL{U^) 

jj{2) ^ ^ 1 ^( 1 ) ^ ‘^AtLiU^^'^) 

4 4 4 

un+i ^ \jjn ^ 2^(2) ^ lAtL{U^^^) . (13) 

3 3 3 

4. Numerical Tests 

SRHD schemes are usually tested by a series of standard benchmark problems. The most 
important of these are the so-called shock tube Riemann problems. A shock tube test in one 



dimension can be described as follows; a one dimensional pipe is divided into two halves by 
a membrane and the thermodynamic states of each half is specihed. When the membrane 
between the halves is removed, depending on the initial conditions, contact discontinuities, 
shocks and rarefaction waves etc. will result. Analytic solutions of the amplitude and position 
of these features have been derived for both classical and relativistic hydrodynamics. For 
analytical solutions of the relativistic hydrodynamics shock tube problems, see Marti & 
Muller (1994); Thompson (1986); Pons et ah (2000). For a list of some standard shock 
tube tests along a review of the performance of some SRHD schemes, see Marti & Muller 
(2003). In addition to the Riemann problems, we have done some other tests that do not 
have analytic solutions. In each case, the results are compared to the literature. 

In the following two sections, we present the results of the benchmark tests in one and 
two dimensions. Before we present our results, we list in Table 1, the input parameters we 
have used in the tests for contact steepening, flattening and the monotonicity preserving 
steps. We have set 7 = 5/3 in all the computations, unless otherwise specihed. 


4.1. One Dimensional tests 


Problem 1: 

The initial states of the right and left halves of the domain are given by pl = 1-0, pn = 
1.0 ,ul = 0.9 (left) and pr = 10.0,pr = 1.0,ur = 0 (right). This initial condition leads 
to a strong reverse shock. In Figure 1, the results are shown for both 3^'^ and 4*^ order 
reconstructions. We have shown the density, velocity and pressure prohles at t=0.4 and 
compared the numerical results to analytic solutions. Excellent agreement between these 
results and the analytic solutions is noted. Also, these results match closely with those ob¬ 
tained by Lucas-Serano et al. (2004) which were obtained using the PPM. As was the case in 
Lucas-Serano et al. (2004), small oscillations noted behind the shocks disappear completely 
as the CFL number is lowered. However, such oscillations remain more pronounced for the 
4 th Qffigr reconstruction. In Table 2, L} errors for the density for both 3^^^ and 4*^ order 
reconstructions are shown. It shows that the order of the scheme is approximately one, as 
is expected for solutions with shocks. 

Problem 2: 

The initial states of the right and left halves of the domain are pl = 10.0, pl = 1.0, ul = 
—0.6 (left) andpR = 10.0, pr = 1.0, ur = 0.5 (right). In Figure 2, the results using the 3'"'^ and 
4 th WENO reconstructions are presented. Once again, we note that the shock positions 
and velocities are captured well by both reconstruction schemes. Unlike the previous test. 
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post shock oscillations are not seen in either of the reconstructions. Direct comparison with 
Lucas-Serano et ah (2004), show good agreements as well. Table 3 presents the errors of 
the density for both 3’’'^ and 4*^ order reconstructions. The order is approximately one for 
both cases except for very high grid spacings. 

Problem 3: 

The initial states of the right and left handed halves of the domain are given by pl = 
13.3, pl = 10.0, ul = 0 (left) and pr = 0, pr = 1.0, ur = 0 (right). We expect a shock wave, 
a contact discontinuity and a rarefaction wave from this test. Figure 3 shows the results of 
the computation. There is good agreement between the computed and analytic solutions. 
Direct comparison with Lucas-Serano et al. (2004) and Zhang & MacFadyen (2005) show 
good agreements as well. In Table 4, the errors of the density for both for 3’’'^ and 4*^ 
order reconstructions are shown. Although irregular, these values show similar trends and 
magnitudes to those obtained by others (see Lucas-Serano et al. (2004), Zhang & MacFadyen 
(2005)). We believe the error could be optimized by appropriately hne tuning the smoothing 
and flattening parameters. 

Problem 4: 

The initial states of the right and left handed halves of the domain are given by pl = 
1000.0, Pl = 1.0, ul = 0 (left) and pr = 0.01, pr = 1.0, ur = 0 (right). The initial condition 
gives rise to a right-moving shock, a left-moving rarefaction wave and a contact discontinuity 
in between. This is a fairly demanding test due to the 10® order pressure ratio. In Figure 4, 
we show the results of our computation. As expected the 3'"'^ order scheme is more dissipative 
than the 4*^ order scheme. Once again, a direct comparison between the results and analytic 
solutions show comparable agreements to the same tests done with other schemes (see Lucas- 
Serano et al. (2004), MacFadyen and Zhang (2005)). The errors (Table 5) of the density 
also show similar results. 

Problem 5 (Shock Reflection Test) We consider now an even more severe test for 
the scheme. The shock reflection test consists of an ultra-cold relativistic wind hitting a 
solid wall and a reflected shock wave propagating backward leaving a static region of hot 
gas. The Lorenz factor of the wind is set by the precision level of the code. Hence, winds 
at very high speeds are considered. To directly compare our results with that of Del Zanna 
& Bucciantini (2002), Lucas-Serano et al. (2004), and Zhang & MacFadyen (2005), we have 
chosen the following initial conditions; p = .01, p = 1.0, n = 0.99999, n = 0.9999999999 for 
the first and second tests, respectively. The solid wall is placed at the right edge of the 
computational domain. The post shock density is an increasing function of the initial flow 
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velocity. The compression ratio a = P 2 /P 1 , satisfies 


a = 


7 + 1 
7-1 


7 


7 


-62 


(14) 


where, 62 = hhi — 1 and Wi is the initial Lorentz factor. We have used 7 = 4/3 for this test. 
The results are shown in Figures 5 and 6 . Consider first Figure 5, which corresponds to a 
Lorenz boost factor, W = 224. Direct comparison between our results and those of Lucas- 
Serano et ah (2004), and Del Zanna & Bucciantini (2002), show good qualitative agreement. 
Note that both 3'"'^ and 4*^ order reconstruction schemes capture the shock wave without 
oscillations, in contrast to the case of Del Zanna and Bucciantini (2002), where the 3’’'^ order 
CENO scheme fails. Next, consider the case of hF = 70710.675 (Figure 6 ). Comparing 
our results to those of Zhang & MacFadyen (2005), there is excellent qualitative agreement. 
These results show that the scheme is robust against both 3'"'^ and 4*^ order reconstructions 
for ultra-relativistic flows. 


Problem 6: Mixed Solution test 

For the final one dimensional test, we followed Del Zanna and Bucciantini (2002), and 
considered a test in which oscillatory and smooth solutions arise close together simultane¬ 
ously. The initial configuration is given by; 


Pl = 50.0, pl = 5.0, Ui = O(left) , (15) 

Pr = 5.0, pr = 2.0 + 0.3 sin 50a;, ur = O(right) . (16) 

The solution consists of the interaction between a blast wave and a density wave. The 
results are shown in Figure 7, which shows the density profile for both the 3'"'^ and 4*^ order 
reconstructions simultaneously. We see that the 4*^ order scheme gives more oscillatory 
results than 3'"'^ order WENO. Direct comparison with Del Zanna & Bucciantini (2002), 
show good agreement. 


5. 2-dimensional Relativistic Hydrodynamics 

The one dimensional scheme is extended to two dimensions in a straightforward manner 
using a dimensionally unsplit scheme for advancing the solution that was described earlier. 
The two dimensional reconstruction is done using dimensional splitting. For the multidimen¬ 
sional tests, we consider the a two dimensional Riemann problem, a symmetric blast wave 
and a relativistic jet. All of these tests are compared to previous results. 
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5.1. Two-dimensional shock tube test 


Following Del Zanna & Bucciantini (2002), Lucas-Serano et al. (2004), and Zhang 
and MacFadyen (2005), consider the following two-dimensional Riemann Problem. A two 
dimensional square region is divided into four quadrants and the initial conditions given by 


{P,Vx,Vy,p)^^ = 
{P,Vx,Vy,p)^^ = 
{P,Vx,Vy,pf^ = 
(p, Va:,Vy,pf^ = 


( 0 . 1 , 0 . 0 , 0 . 0 , 0 . 01 ) , 
(0.1,0.99,0.0,1.0) , 
(0.5, 0.0, 0.0,1.0) , 
(0.1,0.0,0.99,1.0) . 


(17) 


Outflow boundary conditions are used throughout the domain and the adiabatic gas index 
is set to 7 = 5/3. The resulting pattern is shown in Figure 8 which were computed using 
the 4*^ order WENO and PPM (from Lucas-Serano et al. (2004)) reconstruction schemes. 
The solutions are shown at time time t = 0.4. The PPM based calculations have been 
presented here for the sake of direct comparison. The salient features of these results are as 
follows. Note that the WENO based scheme gives sharper contact discontinuities than the 
PPM scheme for the same grid resolution. Around the head of the shock, we have been able 
to resolve both the bow shock and the inner shock structure. Around the tails, where the 
shock merge, we are able resolve the jet-like structure seen due to the merging of the two 
discontinuities. 


5.2. Spherically symmetric blast-wave test 

In this test, a cylindrically symmetric blast wave is studied. The set-up of the problem 
is as follows. The domain is [—1,1] X [-1, 1] with outflow boundary conditions. The initial 
conhguration is given by 

(p,na;Wy,p) = (l-,0.,0., 1000.) if < (.4)^ 

(p, Uj;, ny,p) = (1., 0., 0., 1.0) Otherwise . (18) 

This problem does not have an analytic solution so the results are compared to a one dimen¬ 
sional simulation. The results are also compared to Del Zanna & Bucciantini (2002), and 
Zhang & MacFadyen (2005). The initial condition given above give rise to a cylindrically 
symmetric blast wave moving outward, a rarefaction wave moving inward and and a contact 
discontinuity. A 4*^ order WENO data reconstruction algorithm was used for this test. In 
Fig. 9, the test results are shown at time t = 0.4. Direct comparison with the above men¬ 
tioned references show excellent agreement. Cylindrical symmetry is also well preserved in 
the computations. 
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5.3. 2D Relativistic Jet 

The final multidimensional test undertaken for our algorithm is that of relativistic jets in 
planar geometry. The dynamics and morphology of jets have been studied by many authors 
(see Zhang & MacFadyen (2005), and references therein). This is a standard benchmark for 
testing any new relativistic HRSC scheme. To test the ability of our codes to capture the 
salient features of relativistic jets, we have considered the following scenario. Our tests are 
done in planar geometry. We consider a rectangular domain of length [0,45] x [—25,25]. 
Through a small opening along the axis of symmetry we inject a relativistic fluid at Vx = 
0.99c, Vy = 0.0. The pressure between the inflow and the ambient medium is considered to 
be at equilibrium. And we take a density ratio of the inflow and ambient medium of 0.01. 
Reflective boundary conditions are used along the axis-of symmetry. Otherwise outflow 
boundary conditions are used elsewhere. 

Figs. 10-11 progression of the jets at times t = 5, t = 25, t = 50 and t = 75 respectively. 
Comparing our results to Lucas-Serano (2004), Del Za un a & Bucciantini (2002), and Zhang 
& MacFadyen (2005), we find good qualitative agreement. The jets capture the expected 
morphology quite well. There is a supersonic jet that extends from the nozzle and impacts 
with the ambient medium, at which point we note a terminal planar shock, a contact discon¬ 
tinuity and a bow shock. There is a cocoon that separates the ambient medium from the jet 
material by a contact discontinuity where we note Kelvin-Helmholtz instabilities. We also 
note the absence of carbunckles, that is purely a numerical artifact. 

Figs. 12-13, show snapshots of the same calculations using a slightly narrower opening 
of the in flowing jets. Compared to the previous case, the jets are more collimated in this 
case. 


6. Summary and Conclusion 

We have developed a new multidimensional relativistic hydrodynamics code based on 
the semidiscrete central and WENO reconstruction approach, and some elements of the PPM 
method. This is the first such multidimensional relativistic code combining these elements. A 
dimensionally unsplit scheme was used to advance the solution. We also carried out a number 
of one and two dimensional tests of the code using both 3’’'^ and 4*^ order reconstruction 
methods. The test results are comparable to a number of other codes currently being used 
in computational astrophysics. In all cases, the tests showed good agreements with codes 
currently in use. 

The work presented here should be considered as the first few steps towards the de- 
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velopment of a multi-purpose relativistic hydrodynamics code. Some of the developments 
mentioned above on semidiscrete central schemes could well extend the performance of our 
code. These include the genuinely multidimensional formulation, the formulation on un¬ 
structured grid or higher order WENO reconstructions. From the point of applications, we 
are working to extend our code to three dimensions, provide it with adaptive capabilities 
and apply it for parallel architecture using MPI. 
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Test 

Ko 

7^(1) 

7^(2) 

6(1) 

o;(i) 

o;(2) 

e(2) 

Test 1 

1.0 

5.0 

0.05 

0.7 

0.52 

10.0 

5. 

Test 2 

1.0 

5.0 

0.05 

0.7 

0.52 

10.0 

5.0 

Test 3 

1.0 

5.0 

0.05 

0.7 

0.52 

10.0 

5.0 

Test 4 

1.0 

5.0 

0.05 

0.7 

0.52 

10.0 

5.0 

Test 5 

1.0 

5.0 

0.1 

0.1 

0.52 

10.0 

0.0001 


Table 1: Values of the flattening, steepening and monotonicity preserving parameters used 
in one dimensional shock tube tests presented in Sec. 5.1 


N 

iyi-error (4*i® order) 

Rate 

Ri-error (3'"'^ order) 

Rate 

40 

2.82E-1 

- 

3.17E-1 

- 

80 

1.33E-1 

1.08 

1.55E-1 

1.06 

160 

6.77E-2 

.9466 

8.04E-2 

0.96 

320 

3.4E-2 

1.02 

4.22E-2 

0.94 

640 

1.75E-2 

0.98 

2.32E-2 

0.92 


Table 2: errors of density in Test 1, Sec. 5.1. Shown are the errors in density for both 3'’'^ 

and 4*^ order reconstruction. 


N 

Ri-error (4*^ order) 

Rate 

Ri-error order) 

Rate 

40 

2.53E-1 

- 

3.76E-1 

- 

80 

1.66E-1 

0.63 

2.51E-1 

0.61 

160 

7.31E-2 

1.22 

1.26E-1 

0.99 

320 

3.77E-2 

0.98 

6.55E-2 

0.95 

640 

1.92E-2 

1.15 

3.33E-2 

0.99 


Table 3: errors for Test 2, Sec. 5.1. Shown are the errors in density for both 3’’'^ and 4*^ 

order reconstruction. 
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N 

L^-error (4*^ order) 

Rate 

L^-error order) 

Rate 

40 

4.40E-1 

- 

4.98E-1 

- 

80 

2.28E-1 

0.95 

2.80E-1 

0.84 

160 

1.30E-1 

0.81 

1.49E-1 

0.92 

320 

6.19E-2 

1.07 

7.23E-2 

1.04 

640 

3.23E-2 

0.94 

3.77E-2 

0.94 


Table 4: errors for Test 3, Sec. 5.1. Shown are the errors in density for both 3’’'^ and 4*^ 

order reconstruction. 


N 

L^-error (4*^ order) 

Rate 

L^-error (3^^^ order) 

Rate 

80 

2.79E-1 

- 

2.81E-1 

- 

160 

1.85E-1 

0.51 

1.85E-1 

0.6 

320 

1.51E-1 

0.29 

1.50E-1 

0.3 

640 

l.OlE-1 

0.57 

9.42E-2 

0.66 


Table 5: errors for Test 3, Sec. 5.1. Shown are the errors in density for both 3’’'^ and 4*^ 

order reconstruction. 
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Fig. 1.— Density, Pressure and velocity profiles for test 1 (Sec. 5.1) at t=0.4 using 3^*^ (top) 
and 4*^ (bottom) order WENO reconstructions and n=400. 
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Fig. 2.— Density, Pressure and velocity profiles for test 2 (Sec. 5.1) at t=0.4 using 3^*^ (top) 
and 4*^ (bottom) order WENO reconstructions and n=400. 
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Fig. 3.— Density, Pressure and velocity profiles for test 3 (Sec. 5.1) at t=0.35 using 3^*^ 
(top) and 4*^ (bottom) order WENO reconstructions and n=400. 
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Fig. 4.— Density, Pressure and velocity profiles for test 4 (Sec. 5.1) at t=0.4 using 3^*^ (top) 
and 4*^ (bottom) order WENO reconstructions and n=400. 
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Fig. 5.— Density, Pressure and velocity profiles for test 5 (Sec. 5.1), v= .99999 using 3^*^ 
(top) and 4*^ (bottom) order WENO reconstructions and n=400. 
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Fig. 6.— Density, Pressure and velocity profiles for test 5 (Sec. 5.1), v= .9999999999 using 
3'”'^ (top) and 4*^ (bottom) order WENO reconstructions and n=400. 
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Fig. 7.— Mixed soulution test (test 5, Sec. 5.1) for 3'"'^ and 4*^ order WENO reconstruction 
at t= .36. Shown ar the density prohles. 
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Fig. 8.— 2D Riemann problem, (Sec. 6.1). Logarithmic density plots are shown using 
our scheme (top) using the 4*^ order WENO schme and that using the PPM reconstruction 
method (bottom) of Lucas-Serano et ah (2004). Solution showm at t=0.4 and n=400 
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Fig. 9.— Cylindrically symmetric blast wave (Sec. 6.2) using 4*^ order WENO. Solution at 
t=0.4 and n=250. Solid lines represent one dimensional calculations. 
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